% T1_mCINE_IR: Estimates parameters and errorbars for T1 recovery in an IR sequence 
% Used signal model: S = A*[ 1-B*exp(-R1*TI) + exp(-R1*TR) ] 
% Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)

% After executing a file browser pops up, where the image to be mapped
% needs to be selected. All images should be in the same folder. The
% selected image should be from the first cardiac cycle.

% Two ROIs need to be drawn: First a ROI with only noise and second the ROI
% 
% where the mapping has to be performed

% There is a switch 'roibased' for performing voxel/roi based mappings at the top of
% this script.

% There is a variable 'repeatphase', which should be the number of images
% after which the same cardiac phase is reached in the CINE acquisition.

% Output are dicom images of the R1map, R1errormap, T1map and T1errormap, saved
% in the directory of the first image of the serie. 
% IMPORTANT: The values of the R2 and R2 errormap are multiplied with 10000 to get into sensible gray values.
% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.

% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot,
% University of Antwerp, 13-8-2007.

clear all

% switch for voxel/roi/roimean type of estimation
% TODO: fix error with option 1
roibased = 2; % 0=voxel, 1=roi all voxels, 2=roi mean

% choose if you want a T1 or a R1 map overlay image
overlaytype = 0; % 0=T1, 1=R1

% fill after how many images the same cardiac phase is reached
repeatphase = 3;

% which should be the first point to take into account
fitfrom = 5; 

% which is approximately the first point that is true positive
lsfitfrom = 12;

% location of mask
maskloc = 'phantommask.mat';

% [file,chosendir]=uigetfile('*.dcm');
[file,chosendir]=uigetfile('*.dcm');
imagenr = str2num(file(1:end-4));
cd(chosendir);
d=dir;
nfol=(length(d)-2)/repeatphase;
cd ..
outputdir = pwd; %save maps in dir of chosen image

%read in the image that is picked by the user
chosenimage = dicomread(dicominfo(fullfile(chosendir,file)));

%get the dicom info
for i=1:nfol   
    %henk orig fitdata(i).folder=char(d((i-1)*repeatphase+2+imagenr).name); 
    fitdata(i).folder=char(d((i-1)*repeatphase+2+1).name); 
    info = dicominfo(fullfile(chosendir,fitdata(i).folder));
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime+99999; %TR is not correctly in header, should look it up. It's long compared to T1 anyway, so almost neglectible
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
    fitdata(i).TI=info.TriggerTime;
end

%load masks
load(maskloc);

%determine sigma of the noise
%henk orig imshow(fitdata(2).image,[90 570]);
% imshow(double(chosenimage),[])
% M=roipoly;
% close;
load('C:\DicomQRPush\t1_mapping_gado_phantom_henk_mcine_1985\noisephantom8.mat')

% read in data to fitdata struct
for i=1:nfol
    D=fitdata(i).image.*M;
    sdyd(i)=std(nonzeros(D));
    maxsd=max(sdyd);
end

clear D;

% conversion factor from standard to rayleigh distribution
CRLBsigma=1.527*maxsd;
% CRLBsigma = 49;

%draw masks
    for k=1:nfol
        TI(k)=double(fitdata(k).TI);
        TR(k)=double(fitdata(k).TR);
    end

    TI=TI';
    TR=TR';
   
    
    if ~roibased
        [row,column]=find(fitdata(2).image);
        nrvoxels=size(row,1);
    end;

    dims=size(fitdata(1).image);
%     T1=zeros(dims(1),dims(2));
%     R1=T1;
%     A=T1;
%     B=T1;
%     STDCRT1=T1;
%     STDCRR1=T1;
%     STDCRB=T1;
%     STDCRA=T1;
    
     
    %set optimization options
    opt = optimset('fminunc');
    opt = optimset(opt,'Diagnostics','off','LargeScale','off','gradObj','on','Display','off','MaxIter',80,'Hessian','off','TolFun',1e-12,'Tolx',1e-8);
    s = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,100],'Startpoint',[700, 2, 0.001],'Maxiter',50,'Display','off','TolFun',10^-5,'TolX',10^-5);
    fiteq = (['a*(1-b*exp(-x*c)+ exp(-(' num2str(TR(1)) '+x)*c))' ]);
    % fiteq = (['(a-b*exp(-x*c))']);
    f = fittype(fiteq,'options',s);
    fitrange = lsfitfrom:size(TI)-(fitfrom-1);
    
for p=1:size(masks,2)
    
    for i=1:nfol
        fitdatamask(i).image=fitdata(i).image.*masks(p).mask;

    end

    %initialize&allocate





    %logistics for roibased all voxels
     if roibased == 1
         yd=double(fitdata(1).image);
         for k=2:nfol
                yd=[yd double(fitdata(k).image)];
         end
         yd=reshape(yd,size(yd,1)*size(yd,2),1);
         TI=repmat(TI,1,size(yd,1)/(nfol));
         TI=TI';
         TI=reshape(TI,size(yd,1),1);
         [nonzerox,nonzeroy]=find(yd);
         yd=yd(nonzerox,1);
         TI=TI(nonzerox,1);
         nrvoxels=1;
         fitfrom = (fitfrom-1)*size(nonzeros(M),1)+1;

    % logistics for roibased mean of voxels
     clear yd;
     yd=zeros(size(TI,1),size(TI,2));
     elseif roibased ==2
        yd=mean(double(nonzeros(fitdatamask(1).image)));
             for k=2:nfol
                yd=[yd mean(double(nonzeros(fitdatamask(k).image)))];
             end  
            yd=yd';
            nrvoxels = 1;      
     end

 % estimation loop
    % voxelbased fit
for i=1:nrvoxels
%     disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels(1))]);
        if ~roibased
            if fitdata(1).image(row(i),column(i)) > 0;
                for k=1:nfol
                    yd(k)=double(fitdata(k).image(row(i),column(i)));
                end
                
                [LS,ML] = T1IRAbsComputePar(yd(fitfrom:end),TI(fitfrom:end),0,CRLBsigma,opt,s,f, TR(fitfrom:end),fitrange);
                [CR] = T1IRAbsCramerRao(ML,TI(fitfrom:end),CRLBsigma, 10000.*TI(fitfrom:end));

                A(row(i),column(i))=ML(1,1);
                B(row(i),column(i)) = ML(2,1);
                T1(row(i),column(i)) = 1/ML(3,1);
                R1(row(i),column(i)) = ML(3,1);

                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                else
                    STDCRA(row(i),column(i)) = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB(row(i),column(i)) = sqrt(CR(2,2));
                else
                    STDCRB(row(i),column(i)) = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1(row(i),column(i)) = sqrt(CR(3,3));
                    STDCRT1(row(i),column(i)) = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1(row(i),column(i)) = 5;
                    STDCRT1(row(i),column(i)) = 10000;
                end
            else
                A(row(i),column(i))= 0;
                B(row(i),column(i)) = 0;
                T1(row(i),column(i)) = 0;
                R1(row(i),column(i)) = 0;

                STDCRR1(row(i),column(i)) = 0;
                STDCRT1(row(i),column(i)) = 0;
            end
        %roibased fit
        else 
                [LS,ML]=T1IRAbsComputePar(yd(fitfrom:end),TI(fitfrom:end),0,CRLBsigma,opt,s,f, TR(fitfrom:end),fitrange);
                [CR]=T1IRAbsCramerRao(ML,TI(fitfrom:end),CRLBsigma, 10000.*TI(fitfrom:end));
                A(p) = ML(1,1);
                B(p) = ML(2,1);
                T1(p) = 1/ML(3,1);
                R1(p) = ML(3,1);
                Als(p) = LS(1,1);
                Bls(p) = LS(2,1);
                T1ls(p) = 1/LS(3,1);
                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(p) = sqrt(CR(1,1));
                else
                    STDCRA(p) = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB(p) = sqrt(CR(2,2));
                else
                    STDCRB(p) = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1(p) = sqrt(CR(3,3));
                    STDCRT1(p) = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1(p) = 5;
                    STDCRT1(p) = 10000;
                end
        end
            
    end
end

% Set fit parameters.
fiteq2 = ('a-b*exp(-x*c)');
s2 = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,30000],'Startpoint',[500, 100, 0],'Maxiter',500,'Display','off','TolFun',10^-15,'TolX',10^-15);
f2 = fittype(fiteq2,'options',s2);
yd2 = yd(9:end);
TI2 = TI(9:end);

% Do the fitting
[c22,gof,output] = fit(TI2,yd2,f2,s2);
% output.exitflag
A2=c22.a ;
B2=c22.b;
R12=c22.c;


%output message
% disp(['Estimation finished. Used signal model: S = A*[ 1-B*exp(-R1*TI) + exp(-R1*TR) ] '  ])
% disp(['Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)'  ])
% disp(['Median A = ' num2str(median(nonzeros(A))) ', median B = ' num2str(median(nonzeros(B))) ', median T1 = ' num2str(median(nonzeros(T1)))  ])

%write result in dcm files. R1 values*10000 to get in a sensible range
if ~roibased
    dicomwrite(int16(10000*R1),[outputdir,filesep,file(1:end-4),'_R1map','.dcm'],fitdata(1).info);
    dicomwrite(int16(T1),[outputdir,filesep,file(1:end-4),'_T1map','.dcm'],fitdata(1).info);
    %dicomwrite(int16(A),[maskdir,'\',maskfile,'Amap','.dcm'],fitdata(1).info);
    %dicomwrite(int16(B),[maskdir,'\',maskfile,'Bmap','.dcm'],fitdata(1).info);
    dicomwrite(int16(10000*STDCRR1),[outputdir,filesep,file(1:end-4),'_R1errormap','.dcm'],fitdata(1).info);
    dicomwrite(int16(STDCRT1),[outputdir,filesep,file(1:end-4),'_T1errormap','.dcm'],fitdata(1).info);
    
    disp(['Resulting maps are stored in ' outputdir ])
    
    if ~overlaytype
        %T1 overlay image. The "[# #]" values are the range of the T1 map.
        figure('Name','Overlay with T1 map (milliseconds)', 'position', [800 600 500 500]);
        sc(T1,[00  1500], jet,sc(chosenimage, [0 max(max(chosenimage))],gray),T1==0);
        colorbar('south','XColor','white');
    else
        %R1 overlay image. The "[# #]" values are the range of the R1 map.
        figure('Name','Overlay with R1 map (R1*10000, 1/milliseconds)', 'position', [800 600 500 500]);
        sc(R1*10000,[0 100], jet,sc(chosenimage, [0 max(max(chosenimage))],gray),R1==0);
        colorbar('south','XColor','white');
    end

end

%display plots if roibased
if roibased
    scatter(TI,yd,'*k')
    hold on;
%     scatter(TI(1:9),yd(1:9),'*c')
    ti=floor(min(TI)):1:round(max(TI));
    yy=abs(ML(1,1)*(1-ML(2,1)*exp(-ti*ML(3,1))+exp(-(TR(1))*ML(3,1))));
    plot(ti,yy,'--r');
%     disp(['Model = Abs(A - B*exp(-R1*TI)'])
%     disp(['Results: value +- standard deviation:'])
%     disp(['R1 = ' num2str(ML(3,1)) ' +- ' num2str(sqrt(CR(3,3)))])
%     disp(['T1* = ' num2str(1/ML(3,1)) ' +- ' num2str(sqrt(CR(3,3))/(ML(3,1)^2))])
%     disp(['A = ' num2str(ML(1,1)) ' +- ' num2str(sqrt(CR(1,1)))])
%     disp(['B = ' num2str(ML(2,1)) ' +- ' num2str(sqrt(CR(2,2)))])
%     disp(['T1 [T1* (B/A - 1)] = ' num2str((1/R12)*(B2/A2 - 1))])
T1
end;

    
